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Abstract 

We present lattice calculations for the ground state energy of dilute neutron matter at 
next-to-leading order in chiral effective field theory. This study follows a series of recent 
papers on low-energy nuclear physics using chiral effective field theory on the lattice. In this 
work we introduce an improved spin- and isospin-projected leading-order action which allows 
for a perturbative treatment of corrections at next-to-leading order and smaller estimated 
errors. Using auxiliary fields and Euclidean-time projection Monte Carlo, we compute the 
ground state of 8, 12, and 16 neutrons in a periodic cube, covering a density range from 2% 
to 10% of normal nuclear density. 
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I. INTRODUCTION 



Chiral effective field theory for low-energy nucleons on the lattice has been investigated 
in several recent papers. In Ref. l| chiral effective field theory was considered at leading 
order (LO) using two different lattice actions. These actions, LOi and LO2, each include 
the leading-order interactions in Weinberg's power counting scheme {2!, 3]. The difference 
is that in LOi the nucleon-nucleon "contact" interactions are point-like while in LO2 they 
are smeared using a Gaussian function. These smeared interactions in LO2 were used to 
better reproduce S-wave phase shifts for nucleon momenta up to the pion mass. If the 
effective field theory expansion is converging properly then low-energy physical observables 
computed using LOi and LO2 should agree up to differences of the size of next-to-leading 
order (NLO) corrections. Similarly when NLO corrections are included, agreement should 
be comparable to corrections at next-to- next-to-leading order (NNLO). 

n n 

In Ref. the spherical wall method p\ was used to calculate nucleon-nucleon scattering 
phase shifts and the S-D mixing angle for LOi and LO2 at spatial lattice spacing a = (100 
MeV)~^ and temporal lattice spacing at = (70 MeV)^^. In a companion paper the 
same LO2 lattice action was reproduced using auxiliary fields, and the ground state energy 
of dilute neutron matter was calculated using projection Monte Carlo at densities ranging 
from 2% to 8% of normal nuclear density. For each Monte Carlo configuration next-to- 
leading order corrections were computed using first-order perturbation theory. Simulations 
using the lattice action LOi were also attempted, however strong complex phase oscillations 
prevented an accurate calculation. 

Ground state energy results using the LO2 action at leading order and next-to-leading 
order are shown in Fig. [H The energy is plotted as a fraction of the ground state energy for 
non- interacting neutrons at the same Fermi momentum kp- For kp less than 100 MeV, the 
difference between results at leading order and next-to-leading order is small enough that the 
convergence of the effective theory appears reliable. However for kp greater than 100 MeV 
the difference is relatively large, and the perturbative treatment of NLO corrections seems 
questionable. The analysis in Ref. jo] found that much of the difference between the LO2 and 
NLO2 results could be ascribed to differences in the P-wave phase shifts. Although helpful 
in S-wave channels, the Gaussian smearing used in LO2 produces unphysical attractive 
forces in each P-wave channel which must be cancelled at next-to-leading order. In this 
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FIG. 1: Ground state energy ratio Eq/Eq^^ for dilute neutron matter versus Fermi momentum kp 
for LO2 and NLO2 



paper we introduce a new leading-order action LO3 that solves this problem. The new 
action equals LO2 in each S'-wave channel but matches LOi in each P-wave channel. We 
construct the new action using projection operators for the spin-singlet /isospin-triplet and 
spin-triplet / isospin-singlet channels. 

The paper is organized as follows. We first review the effective potential for chiral 
effective field theory up to next-to-leading order and simplifications that can be made at low 
cutoff momentum. We also summarize the lattice transfer matrix formalism for LOi and 
LO2. The new action LO3 is then introduced and phase shifts and the S-D mixing angle 
are computed up to next-to-leading order. After this we rewrite the LO3 transfer matrix in 
terms of one-body interactions with auxiliary fields. This allows us to simulate the ground 
state of the many-neutron system up to next-to- leading order using projection Monte Carlo. 
We compare the new results obtained using LO3 and NLO3 with the LO2 and NLO2 results 



from Ref. 



6j] and other published data in the literature. We also analyze the ground state 



energy ratio Eq/Eq as an expansion near the unitarity limit. 
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II. CHIRAL EFFECTIVE FIELD THEORY 



A. Effective potential 



In our notation q denotes the t-channel momentum transfer for nucleon-nucleon scatter- 
ing and k is the w-channel exchanged momentum transfer. At leading order in the Weinberg 
power-counting scheme the nucleon-nucleon effective potential includes two indepen- 

dent contact terms and instantaneous one-pion exchange (OPEP), 



1/(0) = Cs + Ct (<?i ■ ^2) 



V 



OPEP 
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g2 _|_ ^2 



(1) 

(2) 
(3) 



The vector arrow in a signifies the three- vector index for spin, and the boldface for r signifies 
the three-vector index for isospin. For physical constants we take m = 938.92 MeV as the 
nucleon mass, = 138.08 MeV as the pion mass, = 93 MeV as the pion decay constant, 
and qa = 1-26 as the nucleon axial charge. 

At next-to-leading order the effective potential has seven independent contact terms car- 
rying two powers of momentum, corrections to the two LO contact terms, and instantaneous 



two-pion exchange (TPEP) 
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NLO 



Following the notation of Ref. 10|, lUl] we have 
Vio + Ar(°) + V(^) + V;^li^. (4) 



The NLO contact interactions are given by 

Ay(°) = ACs + ACt (ai ■ 0^2) 

1^(2) = Cig2 + C2k^ + (C3g2 + c^p^ (ai • ^2) + ^C^^ + ^2) -{qxk 
+ Cq ((Tl ■ q) {^2 ■ q) + (^ai ■ Pj (a2 ■ Pj , 

and the NLO two-pion exchange potential is 

HQ 



T/TPEP _ _Jj_^J]jLT (n\ 

^^LO " 3847r2/4^^^^ 

TV 



Am 



I (5,1 - Ag\ -l)+q^ (23,1 - 10.1 - l) + ^^2 



^.L{q) [{q-Sr){q-a2)-q'{a,.a2)] 



(5) 



(6) 



(7) 
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where 



Recent reviews of chiral effective field theory can be found in Ref. [14 



(8) 



15, 



ll7|. 



B. Power counting and cutoff momentum 

There have been a number of studies on the short-distance behavior of the one-pion 
exchange potential and the consistency of the Weinberg power-counting scheme. An alter- 



native scheme known as KSW power counting was proposed [18|, ll9|, |20|]. This scheme is 
based on a perturbative treatment of the one-pion exchange potential and allows for system- 
atic control of ultraviolet divergences in the effective theory. Unfortunately convergence at 
higher orders was found to be poor in some partial waves for momenta comparable to the 



pion mass 



2l|. 



Other power counting alternatives have also been proposed. In one scheme the lead- 
ing /~^r~'^ short-distance singularity is treated non-perturbatively while the remainder of 



the one-pion potential is introduced as a perturbative expansion in powers of [22|. 
More recently a different power counting modification was proposed in which one-pion 
exchange is treated non-perturbatively in lower angular momentum channels along with 



higher-derivative counterterms promoted to leading order [23(]. Further investigations of 
this approach in high er p artial waves and power counting with one-pion exchange were 
considered in Ref. 2J, |25|] . 



On the lattice the ultraviolet momentum cutoff is inversely proportional to the lattice 
spacing, A = vr/a. For simple calculations of two-nucleon scattering on the lattice we could 
take any lattice spacing satisfying A ^ vtiT^. However in few- and many-nucleon calculations 
where we use Euclidean-time projection and auxiliary-field Monte Carlo methods, severe 
numerical problems appear when A is very large. In some attractive channels the problem 
is due to spurious deeply-bound states which may appear at sufficiently large A. In other 
channels one faces the problem due to short-range hard-core repulsion at very large A. This 
is manifested as sign or complex phase oscillations which scale exponentially with system 
size and strength of the repulsive interaction. 

To avoid these problems we consider lattice simulations where the cutoff momentum is 
only a few times the pion mass. In this study we take A = 314 MeV ~ 2.3mjr, corresponding 



with a ^ = 100 MeV. For this low cutoff scale the advantages of the alternative power 



counting schemes discussed above are numerically insignificant [26|, and so we use standard 
Weinberg power counting. For nearly all |g| < A we can expand the two-pion exchange 
potential in powers of g^/(4m^), 

+ (10) 



Ami ^ (P' 3 4m^ 

TPEP ''"l ■ ''"2 



(8^1 - 4,1 - 1) + (34,1 - 17,1 - 2) + O (ife) 
[(g ■ o,) (g ■ a^) - (<?! • a^)] [l + O (^)] . (11) 



647r2/.^ 

This expansion fails to converge only for values of q near the cutoff scale A 2.3m^, where 
the effective theory is already problematic due to large cutoff effects of size O icp j There 

is no reason to keep the full non-local structure of V^-^q^ at this lattice spacing. Instead 
we simply use 

Flo = + V^^''''^ (12) 

l^NLo = ^Lo + A\/(°) + F(2), (13) 

where the terms in Eq. ffTTl) with up to two powers of q are absorbed as a redefinition of the 
coefficients AF(°) and F^^), This same approach can be applied to the two-pion exchange 
potential at next-to-next-to-leading order and higher-order n-pion exchange potentials. 

III. LATTICE FORMALISM 
A. Lattice notation 

In this paper we assume exact isospin symmetry and neglect electromagnetic interactions. 
We use n to represent integer-valued lattice vectors on a three-dimensional spatial lattice 
and either p, g*, or k to represent integer- valued momentum lattice vectors. / = 1, 2, 3 
are unit lattice vectors in the spatial directions, a is the spatial lattice spacing, and L is 
the length of the cubic spatial lattice in each direction. We use the Euclidean transfer 
matrix formalism defined in pL] with lattice time step at, and the integer nt labels the time 
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steps. We define at as the ratio between lattice spacings, at = at/ a. Throughout we use 
dimensionless parameters and operators, which correspond with physical values multiplied 
by the appropriate power of a. Final results_^are presented in physical units with the 
corresponding unit stated explicitly. As in [1] the spatial lattice spacing is a = (100 
MeV)^"*^ and temporal lattice spacing is at = (70 MeV)^^. 

We use a and to denote annihilation and creation operators. To avoid confusion we 
make explicit in our lattice notation all spin and isospin indices using 

^0,0 = flT,P' ^0,1 = ^T,"' (l^) 
'^1,0 = O-Up, 0,1,1 = o-i^n- (15) 

The first subscript is for spin and the second subscript is for isospin. We use tj with 
/ = 1,2,3 to represent Pauli matrices acting in isospin space and as with S = 1,2,3 to 
represent Pauli matrices acting in spin space. We also use the letters S and / to denote the 
total spin and total isospin for the two-nucleon system. The intended meaning in each case 
should be clear from the context. We use the eight vertices of a unit cube on the lattice 
to define spatial derivatives. For each spatial direction / = 1,2,3 and any lattice function 
f{n), let 

^ifi^) = \ Yl i-^Y'^'fin + i^), z7=z/ii + /^22 + z/33. (16) 

1^1,1^2,1^3=0,1 

We also define the double spatial derivative along direction I, 

/(n) = f{n + i) + f{n - I) - 2f{n). (17) 



B. Densities and current densities 



We define the local density of nucleons at lattice site n. 



i,i=0,l 



(18) 



This is invariant under Wigner's SU(4) symmetry transforming all spin and isospin degrees 
of freedom [27|. Similarly we define a local spin density for S = 1, 2, 3, 



i,j,i'=0,l 



n] 



(19) 
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isospin density for / = 1,2,3, 

pf'''{n)= Yl «l(^)N,ya..'(n), (20) 

and spin-isospin density for 5", / = 1, 2, 3, 

Psji.n)= ciiji.n)[as\ii,[Ti]..,ai,^j,{n). (21) 

i j'=0,l 

For each static density we also have an associated current density. Similar to the defini- 
tion of the lattice derivative A; in Eq. f|T6|) . we use the eight vertices of a unit cube, 

i7= i/ii + Z/22 + Z/33, (22) 

for i/i, V21 ^3 = 0, 1. Let for / = 1, 2, 3 be the result of reflecting the Z*'^- component of 

V about the center of the cube, 

i7(-/) = z7+(l-2z/;)/. (23) 

Omitting factors of i and 1/m, we can write the /^'^-component of the SU(4)-invariant current 
density as 

nf '"(n) = \ E ^-^y^HA^ + + y)- (24) 

1'1,'^2,!^3=0,1 i,j=0,l 

Similarly the /*'^-component of spin current density is 

</(^) = i E E ^-^r^"<A^ + ^(-0) K].. a,,{n + z7), (25) 
/*^-component of isospin current density is 

</(^) = i E E (-l)'^'+y,(^+^(-0)Kya,,.(n + z7), (26) 

!^lil'2,i^3=0,l j,J,j'=0,l 

and /^'^-component of spin-isospin current density is 

<5}(^) = ^ E E {-ir^'al{n + V{-l))[aslAri],r<^^^ (27) 

ui, 1^2,1/3=0,1 i,j,i',j'=0,l 
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IV. LATTICE ACTIONS 



A. Instantaneous free pion action 

The lattice action for free pions with purely instantaneous propagation is 

Stt7,{t^i) = Q!t(^ + 3) 7r/(n, nt)7r/(n, Ut) - ^ 7r/(n, nt)T:i{n + 1, rit), (28) 

n,nt,I n,nt,I,l 

where tt/ is the pion field labelled with isospin index /. We note that pion fields at different 
time steps rit and n[ are decoupled due to the omission of time derivatives. This decoupling 
among different time steps generates instantaneous propagation in one-pion exchange dia- 
grams and eliminates radiative pion loops. It is convenient to define a rescaled pion field, 

Ti'j{n,nt) = ^TTiifi.nt), (29) 

q^^at{ml + &). (30) 

Then 

S^t,{tt'i) -^]-^ Ti\{n,nt)T:'i{n,nt) -— ^ n'j{n,nt)7r'j{n + i,nt). (31) 

n,nt,I n,nt,I,l 



In momentum space the action is 



SnM) = j^^7^'i{-k,nt)7T'j{k,nt) 



L3 

I,k 



(32) 



The instantaneous pion correlation function at spatial separation n is 

/ J Dn'j exp [-S^^\ 



lj]e-¥^^-"D.(^), (33) 



L3 

k 

where 

DAk) = TTITT- (34) 

l-^E.cos(^) ^ ^ 

B. Transfer matrices for LOi and LO2 

Roughly speaking, the Euclidean-time transfer matrix is the exponential of the Hamil- 
tonian, exp{—HAt), where At equals one temporal lattice spacing. The normal-ordered 
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transfer matrix for non-interacting nucleons is 



Mfree =: exp (-iJfreeOit) : , 



(35) 



where the :: symbols indicate normal ordering. We use the 0(a'^)-improved free lattice 
Hamiltonian, 

49 



free 



12m 



n ij=0,l 



n jj=0,l Z=l,2,3 



+ ^Y.Y. [aUn)a^,An + 2i) + al{n)a,j{n-2i) 

n ij=0,l i=l,2,3 
n i,j=0,lZ=l,2,3 

Let us define the two-derivative pion correlator, 

GsiS2{n) = (^Asyiifi, nt)Asyi{0, n*)^ (no sum on I) 

= ^ E E (-l)'^^M-l)^^^(<(^ + ^-^',^t)7r;(0,< 

With interactions included, the lattice transfer matrix LOi has the form 



(36) 



(37) 



Mloi = : exp i -H,,^at - \cat E [p"''"(n)] ' - \Cj.at E E 

V n In 



-.2 ^,2 



+ E E^^^^^(^i-^^)4:"(^i)p&(^2) 

°J ■K'dir a a T - - 



(38) 



where C is the coefficient of the Wigner SU(4)-invariant contact interaction and C72 is 
the coefficient of the isospin-dependent contact interaction. For the S'-wave there are 
two independent channels corresponding with the spin-singlet /isospin-triplet and the spin- 
triplet/isospin-singlet. To reproduce the physical scattering lengths in each channel we 
set Cs=o,i=i = -5.021 x 10"^ MeV-^ and Cs=i,i=o = -5.714 x 10"^ MeV'^ and use the 
relations 

C = (3C5=o,7=i + Cs=i,i=o) /4, (39) 
C12 = {Cs=o,i=i - Cs=i,i=o) /4. (40) 
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M, 



The LO2 transfer matrix is [l| 

at 



LO2 



: exp ^ -HfreeCtt - 



2L3 



+ 



n2 



(41) 



Si,S2,I ni,n2 

where the momentum-dependent coefficient function /(g) is defined as 



fio) = fo ^exp 



-"Ed 



cos qi ) 



(42) 



and the normahzation factor /o is determined by the condition 



fo 



-y 



exp 



(43) 



-6^ (1 - cosg;, 

As in Ref. [l| we use the value b = 0.6. This gives approximately the correct average effective 
range for the two S-wave channels when C and Cp are tuned to produce the physical S-wave 
scattering lengths. We set ^5=0,/=! = -3.414 x 10"^ MeV^^ and ^5=1,7=0 = -4.780 x 10"^ 
MeV~^ and use the same relations Eq. (l39l) and ( HOl) . The momentum-dependent function 
f{q) produces the Gaussian-smeared "contact" interactions discussed in the introduction. 

The replacement of pointlike interactions in LOi with Gaussian-smeared interactions 
in LO2 is similar to the lattice improvement program of Symanzik used in lattice QCD 



actions 
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29|. There is a conceptual difference however since we are dealing with an 
effective field theory rather than a renormalizable field theory. The higher-order operators 
we consider do not only cancel lattice artifacts but also include higher-order interactions of 
the effective theory. In our lattice calculations the improved leading-order action is treated 
non-perturbatively while higher-order interactions are included as a perturbative expansion. 
The choice of improved action sets a dividing line between perturbative and non-perturbative 
interactions. This dividing line should be immaterial so long as the perturbative expansion 
converges. At any given order, lattice calculations using different improved actions should 
agree up to corrections the size of terms at next order. 



C. Transfer matrix for LO3 

The Gaussian smearing used in LO2 is useful in S'-wave channels but produces unphysical 
attractive forces in P-wave channels. To avoid this problem we introduce a new leading- 
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order action LO3 that equals LO2 in each S'-wave channel but matches LOi in each P- 
wave channel. We multiply the Gaussian-smeared "contact" interactions with projection 
operators for the spin-singlet /isospin-triplet channel, Ps=oj=i, and the spin-triplet /isospin- 
singlet channel, Pg^i j^Q. If we assign labels to the two nucleons, A and B, these projection 
operators are 

Ps..,.. = (i 4 E "l-l) ( J + J E -/-f ) . (44) 

(i + iE-M)(i-iE-^f)- (4S) 

We can define corresponding momentum-dependent density correlations, 

+ ^ ■ J2pf"(^pf"(-^ ■ -4 ■ Y.Psf('i)Psf(-^ ^ (46) 



S=l,/=0 



32 v./r, V 32 

/ SJ 



s 



32 1 V 32 

/ SJ 
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We use Vs=Qj=i and Vs=ij=q to write the leading-order transfer matrix for LO 

MLO3 =: exp ^ -Hi.c^at - j^'^fio) [Cs=o,i=iVs=o,i=i{q) + Cs=i,i=oVs=i,i=o{q)] 

q 

+ E E ^^^^^(^1 - ^^)pti(^^)pti(^2) \ : . (48) 

5i,52,/ ni,n2 J 

The momentum-dependent coefficient function /(g) is the same as defined in Eg. |42] and H3l 



D. Lattice interactions at next-to-leading-order 

n 

The lattice interactions at next-to-leading order were discussed in Ref. We follow 

the same formalism here. We start with the corrections to the leading-order "contact" 

interactions. These NLO interactions are chosen to be point-like rather than smeared 

fl 

operators, and we write the interactions in the same manner as in Ref. 1^, 

A\/ = iAC:5^p'^^''^(n)p'^^''^(n):, (49) 
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= \aCp : ^ pf ''^(n)pf ''^(n) : . (50) 

n,I 

At next-to-leading order there are seven independent contact interactions with two deriva- 
tives. The basis we choose is 

= ■■ Y^P^'^^nP^P^'^^in) :, (51) 

n,l 

Vp^g2 = ~^Cp^g2 : ^pf '''(n)vfpf '"(n) :, (52) 

n,I,l 

= -'^Cs2,g2 : ^pp^inp^Ps'^in) :, (53) 

n,S,l 

Ks2,/2,g2 = -^Cs2j2^g2 : Psf{n)vfpsf{n) :, (54) 

n,S,I,l 

^(q-sr = \C(^g.s)2 : ^sPs'^in) Y2 ^S'pt'^in) :, (55) 

n S S' 

Vp,(q-sr = ^C'/^(g-5)2 : Y2 ^sPsfin) ^ As'Pp^in) :, (56) 

n,7 S S' 

%x5).fe = -^C(.,x5).fc : 5^ [nf ''^(n)A,p";''^(n) + </(n)A,p'^^''^(n)] : . (57) 

n,l,S,l' 

These operators are different from those shown in Eq. ([6]) and allow for a simple projection 
onto different isospin channels. This will be useful later when restricting to the interactions 
of neutrons. 

The V(igxs)-k term corresponds with the continuum interaction 

Ciigxsyk{iqy< {ff^ + ff'')) -k, (58) 

which vanishes unless the total spin is S = 1. The continuum limit of this interaction is 
antisymmetric under the exchange of q and k and is nonzero only for odd parity channels. 
However the lattice interaction V(igx5)-fc does not share this exact t-u channel antisjTiimetry 
at nonzero lattice spacing. Therefore V(^igxs)-k has small lattice artifacts for 5 = 1 in even 
parity channels. We remove this defect by including an exphcit projection onto total isospin 
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/ = 1, 



V, 



1=1 



1=1 



(iqxS)-k C2 {iqxS)-k \ ^ 



n,LS,l' 



^l,S,l' 



n 



(59) 



n,l,S,l',I 

This projection completely eliminates lattice artifacts in the 5* = 1 even parity channels. 



V. SCATTERING RESULTS FOR LO3 AND NLO3 

We measure phase shifts and mixing angles using the spherical wall method 5!]. This 
consists of imposing a hard spherical wall boundary on the relative separation between 
the two nucleons at some chosen radius -Rwaii- Scattering phase shifts are determined 
from the energies of the spherical standing waves, and mixing angles are extracted from 
projections on to spherical harmonics. At next-to-leading order there are nine unknown 
operator coefficients: AC, AC/2, Cq2, Cj2^g2, C52 ^2, C52 j-2 ^2, C(q.5)2, C/2 (^.5)2, and Cj^~^s)-k- 
These nine operator coefficients are fit in the same manner as described in Ref. For 
-Rwaii = 10 + e lattice units, where e is an infinitesimal positive number, we compute energy 
levels for the eight spherical wall modes listed in Table [B The labelling of these modes is 
discussed in Ref. jj]. In addition to these we also consider Q^, the quadrupole moment of 
the deuteron. Qd is a measure of S-D partial wave mixing at low energies and is somewhat 
easier to compute on the lattice than the S-D mixing angle. For each of the nine observables 
we compute derivatives with respect to the nine NLO coefficient operators. By inverting 
the resulting 9x9 Jacobian matrix, we find values for the NLO coefficients needed to match 
each of the nine target values using first-order perturbation theory. The results for the 
operator coefficients are shown in Table [Tll 

With the NLO3 coefficients in hand, we can now calculate lattice phase shifts and mixing 
angles up to next-to-leading order using the spherical wall method. We consider spherical 
walls with radii -R„aii = 10+e, 9 + e, and 8+e lattice units. In order of increasing momentum, 
the lattice data correspond with the first radial excitation for -Rwaii = 10 + e, 9 + e, and 8 + e; 
second radial excitation of -Rwaii = 10 + e, 9 + e, and 8 + e; and so on. The S'-wave phase 
shifts for LO3 and NLO3 versus center-of-mass momentum pcM are shown in Fig. [2l The 
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TABLE I: Results for LO3 and the physical target values 



Splicriccil Wctv6 


Ktpp niirlpons 




PWA93 


l^S'n fMeV) 


0.928 


0.418 


0.407 


3^Sr, (MeV) 


8.535 


6.843 


6.815 


l^S(D)^ (MeV) 


0.928 


-2.225 


-2.225 


3^S(D)^ (MeV) 


8.535 


5.430 


5.675 


2^ Pi (MeV) 


5.691 


5.755 


5.782 


2^P{F)o (MeV) 


5.691 


5.569 


5.584 


2^P{F)i (MeV) 


5.691 


5.754 


5.753 


2^P(F)2 (MeV) 


5.691 


5.684 


5.669 


Qd (fm^) 


N/A 


0.276 


0.286 



TABLE II: Results for NLO3 operator coefficients 



Coefficient 


NLO3 


AC (MeV-2) 


-1.02 X 10"^ 


ACj2 (MeV-2) 


1.03 X 10"^ 


Cq2 (MeV-4) 


2.39 X 10-1° 


Ci2^g2 (MeV-^) 


-4.80 X 10-11 


C52,,2 (MeV-4) 


1.67 X 10-1° 


Cs2j2^q2 (MeV-^) 


-1.03 X 10-1° 


C(,.s)2 (MeV-4) 


-1.43 X 10-1° 


C,2,(,.5)2 (MeV-4) 


1.80 X 10-1° 


Ci=lsyk (MeV-^) 


1.60 X 10-1° 



NLO3 results are in good agreement with partial wave results from Ref. J30(]. 

We plot the S-D mixing parameter Si in the Stapp parameterization [31] in Fig. [31 The 
pairs of points connected by dotted lines indicate pairs of coupled solutions in the spherical 
wall formalism. While there are some deviations from the partial wave data from Ref. 30| , 



the discrepancy is consistent with effects produced by higher-order interactions. As expected 
the 5*- wave results for LO3 are identical with LO2 results in Ref. In fact they agree in all 
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FIG. 2: S'-wave phase shifts versus center-of-mass momentum for LO3 and NLO3. 
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FIG. 3: £1 mixing angle versus center-of-mass momentum for LO3 and NLO3. 

even-L partial wave channels. Results for NLO3 and NLO2 are also close though not exactly 
the same. There are very small differences between the two due to NLO interactions which 
are not completely separable into S'-wave and P-wave terms at nonzero lattice spacing. 

The P-wave phase shifts are shown in Fig. HI We see that the NLO3 results match the 
partial wave data quite accurately. Just as LO3 and LO2 agree in all even-L partial wave 
channels, LO3 and LOi agree in all odd-L partial wave channels. Results for NLO3 and 
NLOi are nearly identical, with only small differences due to the numerical fitting of NLO 
coefficients on the lattice. 
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FIG. 4: P-wave phase shifts versus center-of-mass momentum for LO3 and NLO3. 

VI. AUXILIARY-FIELD FORMALISM FOR NEUTRON MATTER 

So far we have been discussing general systems of low-energy nucleons with both protons 
and neutrons. For computational efficiency we now specialize to the case where all nucleons 
are neutrons. In this case all nucleon-nucleon interactions are in the isospin-triplet channel. 
In the leading-order transfer matrix Mlos '^^^ drop the spin-triplet /isospin-singlet term 
involving Vs=ij=o{(l) and make the simplifying replacements, 

Vs=oj=i{q) - I : P°-'^(g1p°-'^(-g1 : : ^ '"^(^Ip^ ''^(-^1 (60) 

^3,82(^1 - n2)ps^}{ni)Ps^}{n2) ^ Y '^^iSal^l - ^2)P5i'''(^l)P52'"(^2)- 

81,82,1 ni,n2 81,82 ni,n2 

(61) 
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These modifications do not affect the interactions between neutrons and yields the simphfied 
transfer matrix, 



Ml 



LO3 



exp <^ -Hf^eeOit ^77^^ > J{q) 



8L3 



+ 



(62) 



Si,52 ni,ra2 

At next-to-leading order the simplified neutron transfer matrix is 



8L3 



P 



at 



AC ■^^^i^^'-^c 



(9-5)2 



-,2 ^,2 



2n 



(63) 



5i,52, ni,n2 



In the following we use these simplified forms for the leading-order and next-to-leading-order 
transfer matrices. 

The transfer matrices in Eq. ( l62l) and (1631) can be rewritten in terms of one-body inter- 
actions with auxiliary fields. The exact equivalence between lattice formalisms with and 
without auxiliary fields is detailed in HjlsS]. We summarize the results here. 

In neutron-neutron scattering only the neutral pion contributes to one-pion exchange. 
We have been writing the rescaled neutral pion field as TTg, but now we drop the subscript 
"3" and simply write vr'. Let M^^'^^tc' , s, ss) be the leading-order auxiliary-field transfer 
matrix at time step rit, 



M("*)(7r',s,S5) =: exp <j -H^^^^at + -f^Y^s7T'in,nt)pf'' 



n] 



n,S 



n,S 



(64) 
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We can write Mlos normalized integral 

[ Du'DsDss e-^- '-^-''M("')(7r', s, Ss) 

Mlo3 = 7 , (65) 

/ Dn'DsDss e~^^^^ -Sas^ 

where S'i"*^ is the piece of the instantaneous pion action in Eq. ( l28l) containing the neutral 
pion field at time step n^, 

St\'n') = \Y,'^\n,my{n,m) - ^Y.n'{n,n,)n\n + i,n,), (66) 



nl 



and is the auxiliary- field action at time step n^, 



Sts^K^^'^^s) = ]^^s{n,nt)f ^{n-n')s{n' ,nt) + ]^ ^ Ss{n,nt)f ^{n - n')ss{n' ,nt), (67) 

n,n' n,n',S 

with 

The NLO interactions require some additional auxiliary fields. Let 

U^-^\e) = J2^pin,n,)p^'''^{n) + £,,(rl, n,)p^'''^(n) + J2^^sp{n,nt)Asp''''%n) 

n n,S n,S 

+ ^ £Asp^,{^: ^f)Agpg!'°(n) + y^£:v2p(n, ni)V rp"''^(n) 

+ IZ^v2p5(^'^t)^?Ps '''(^) + Xl^n,(n, nt)nf '"(n) + J^£n,,5(^) ^Onf//(n). 

(69) 



ra,Z,5 n,l n,LS 



With these extra fields and linear functional U^^''\e) we define 



M("')(7r',s,S5,e) =: exp \ -H,,,.at + ^|^^A57r'(n, n,)p^'''^(n) 



+ '-^-Cs=,j=,a,Y,^s{n.n,)pi^''{n) + v/^f/("')(^) ^ : • (70) 



We also define the normalized integral, 



/ Dti'DsDss e-^--'-^-*'M("*)(7r', s, ss, e) 
M^'^'\e) = ^ . (71) 



Dtt'DsDss e~'^^'^*'~'^^=*' 
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When all e fields are set to zero we recover Mlos, 

M("')(0) = Mlo3. 



(72) 



To first order the NLO interactions in Mnlos '^^^ t)e written as a sum of bilinear derivatives 
of M("*)(£:) with respect to the e fields at e = 0, 

2 6ep{n,nt) 6ep{n,nt) 

2 5ep{n,nt) 5ey2p{n,nt) 



e=0 



+ ■■•. (73) 

£=0 



VII. EUCLIDEAN-TIME PROJECTION MONTE CARLO 

We extract the properties of the ground state using Euclidean-time projection. We briefly 
summarize the calculation in continuous-time notation before describing the transfer matrix 
calculation at nonzero temporal lattice spacing. Let be a Slater determinant of free- 

particle standing waves in a periodic cube for neutrons. Let -f^LOs be the Hamiltonian 
at leading order and H^lo^ be the Hamiltonian at next-to-leading order. Let -f^su{2)y be 
the same as i^LOs; but with one-pion exchange turned off by setting to zero. As the 
notation suggests, -ffsu(2)^ is invariant under an exact SU(2) intrinsic-spin symmetry. We 
deflne a trial wavefunction 

\^it')) = exp (-i/su(2)/) l*'^^'^) . (74) 

The operator exp (— -ffsu{2)y^') acts as an approximate low-energy filter. In the auxiliary- 
field Monte Carlo calculation this part of the Euclidean-time propagation is positive definite 



for any even number of neutrons invariant under the SU(2) intrinsic-spin symmetry 34 



35, 



36|. With this trial wavefunction we define the amplitude, 

Z{t) = (^(t')l exp {-H^o,t) |vl>(t')) , (75) 
as well as the transient energy at Euclidean time t, 

i?LO.(t) = -|[lnZ(t)]. (76) 
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In the limit of large t, 

lim Ej^osit) = £^o,L03, (77) 

t— >oo 

where -Eq^lOs energy of the lowest eigenstate l^'o) of i^LOs with nonzero inner product 

with 

To compute the expectation value of some general operator O we define 

Zo{t) = exp {-H^Ost/2) O exp {-H^o,t/2) \^{t')) . (78) 

The expectation value of O for |^'o) is given by the large t limit, 

/n^^-(^o|0|*o). (79) 

Corrections to the energy at next-to-leading order can be computed using O — i^NLOa ~Hi^Oz- 
In that case 

J™ = -^0,NLO3 - -E'cLOs) (80) 

t-»oo Z\t) 

where £'o,nlo3 is the ground state energy at next-to-leading order. 
On the lattice we construct \^{t')) using 

|*(0) = (Msu(2)y)^*°|*'^^''>, (81) 

where t' = Lf^at and Lt„ is the number of "outer" time steps. The amplitude Z{t) is defined 
as 

Z(t) = (*(0l(MLO3)^M*(0), (82) 
where t = L^^at and L^^ is the number of "inner" time steps. The transient energy 

Ej^Osit + at/2) (83) 

is given by the ratio of the amplitudes for t and t + at, 

p-BL03(*+«t/2)-at _ + Q^*) (OA) 

~ z{t) ■ ^ ^ 

The ground state energy i?o,L03 equals the asymptotic limit, 

£;o,L03 = lim £^L03 {t + at/2) . (85) 

t— >oo 

We calculate these Euclidean-time projection amphtudes using auxihary fields. For a 
given configuration of auxiliary and pion fields, the contribution to the amplitude Z{t) is 
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proportional to the determinant of an x matrix of one-body amplitudes, where is the 
number of neutrons. Integrations over auxiliary and pion field configurations are computed 



33 



37|. 



using hybrid Monte Carlo. Details of the method can be found in Ref. 

For the ground state energy at next-to-leading order we compute expectation values of 
^NLOa and Mlo3 inserted in the middle of a string of Mloj transfer matrices, 

ZM^.o,it) = (^(i')l (Mlo3)^'^/'Mnlo3 (Mlos)^*'/' , (86) 

ZM,o,it) = m')\ (Mlo3)^*^/'Mlo, (Mlo3)^*^/' \m) ■ (87) 
From the ratio of amplitudes, 

^^^^ = 1 - AENL03 (t)«* + ■ ■ ■ , (88) 

we define the transient NLO energy correction A£'nlo3(^)- The ellipsis denotes terms which 
are beyond first order in the NLO coefficients. The NLO ground state energy -Eo,nlo3 is 
calculated using 

-Eo,NL03 = -So,L03 + lim Ai?o,NL03(^)- (89) 

t — ^oo 

VIII. PRECISION TESTS OF MONTE CARLO SIMULATIONS 

We use the two-neutron system to test the auxiliary-field Monte Carlo simulations. The 
same observables are calculated using both auxiliary-field Monte Carlo and the exact transfer 
matrix without auxiliary fields. We choose a small system so that stochastic errors are small 
enough to expose disagreement at the 0.1% — 1% level. We choose the spatial length of 
lattice to be L = 4 and set the outer time steps Lt^ = 2 and inner time steps Lt- = 2. With 
16 processors we generate a total of about 8 x 10^ hybrid Monte Carlo trajectories. Each 
processor runs completely independent trajectories, and we compute averages and stochastic 
errors by comparing the results of all processors. 

For the first test we choose |\Ef^''^''^ to be a spin-singlet state built from the Slater deter- 
minant of standing waves and \ip2) with 

(0| aij{n) l^i) oc SiflSj^i, (0| aij{n) l^/'s) oc (90) 

For the second test we choose a spin-triplet state with standing waves 

(0| a,,j(n) liPi) oc 6,,ohi cos (0| aij{n) l^s) oc 5i,o5j, i sin (91) 
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TABLE III: Monte Carlo results versus exact transfer matrix calculations for the two-neutron spin 
singlet 5 = and spin triplet 5 = 1. 





5 = (MC) 


5 = (exact) 


5 = 1 (MC) 


5 = 1 (exact) 


Ft ^ (f -i- rv^ /9'\ nVTpVl 
-f^LOs -|- a*/ .ij [ivievj 




— 9 Q1 1 9 


98 '^l'9'> 




Xtc'-f^ [10^ MeV^] 


4.751(5) 


4.7487 


0.0003(7) 







1.580(2) 


1.5789 


-1.025(4) 


-1.0264 




-4.741(6) 


-4.7366 


-1.023(6) 


-1.0264 




-5.788(8) 


-5.7818 


2.51(2) 


2.5533 




0.018(13) 





3.27(16) 


3.4534 


A^NL03(i) [MeV] 


-0.01655(8) 


-0.016440 


-0.2455(9) 


-0.24559 



Comparisons between Monte Carlo results (MC) and exact transfer matrix calculations 
(exact) are shown in Table IIIII The numbers in parentheses are the estimated stochastic 
errors. The agreement between Monte Carlo results and exact transfer calculations is 
consistent with the estimated stochastic errors. 

IX. RESULTS 

We simulate the ground state for = 8, 12, 16 neutrons on periodic cube lattices. For 
= 8 we consider cube lengths L = 4, 5, 6, 7 lattice units. For = 12, we use L = 5,6, 7, 
and for = 16 we use L = Q,7. For each value of A^ and L we fix at either 8 or 10 and 
vary Lt- from 2 up to 12. For |\1>^''°'^^ we take the Slater determinant formed by standing 
waves 

(0| aij{n) |V'2fe+i) oc 6ifi6j^ifk{n), (0| aij{n) |V^2fc+2) oc SiA6j^ifk{n), (92) 

where 

/o(n) = l, /i(n) = cos^, f,{n) =sm'-^, f,{n)=cos'-f^, (93) 

/^(^) = sin2z^, /5(^)=cos^, /6(n) = sin^, /, (n) = cos ^^^^ . (94) 

For A^ = 8 the values of k span the range < A; < 3. For A^ = 12, < A; < 5, and 
for N = 16, < k < 7. For each value of Lf. a total of about 5 x 10^ hybrid Monte 
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Carlo trajectories are generated by 2048 processors, each running completely independent 
trajectories. Averages and stochastic errors are computed by comparing the results of all 
processors. 

Let Eq"^^ be the energy of the ground state for noninteracting neutrons. In Fig. O we 
show the dimensionless ratios 

^free ' ^free ' ^free ' V / 

versus Euclidean time t. These are labelled using the shorthand LO3, ANLO3, and NLO3 
respectively. In addition to the Monte Carlo data we plot the asymptotic expressions, 

A-gNLOa (^) ^ -^0,NLO3 " -^O.LOg j^^-5E t/2 /gy\ 
iTifree ~ rplvee \ / 

EloM + AE^LO,{t) ^ ^0.NLO3 ^^^SE-t Q^SE-tn _ .gg) 
^iree j^tree ^ ' 

The unknown coefficients A and B, energy gap 6E, and ground state energies -E'o.los and 
-£'0,NLO3 are determined by least squares fitting. The e~^^'^ dependence in Eq. fl96|) comes 
from the contribution of low-energy excitations with energy gap 6E above the ground state. 
The e"*^^ */^ dependence in Eq. fp7|) comes from the matrix element of Mnlos between the 
ground state and low-energy excitations at energy gap 6E. 

The resuhs of the asymptotic fits for Eo,loJEI''' and E^, ^1^03/ Eq'^*^ are shown in Table [IVl 
On average the P^r degree of freedom for the fits is around 1. The error estimates for 
£'o,L03 / Eq'^'^ and -Eq.nlos / Eq'^'^ are calculated by explicit simulation. We introduce Gaussian- 
random noise scaled by the error bars of each data point. The fit is repeated many times 
with the random noise included to estimate the one standard-deviation spread in the fit 
parameters. In Table HVl the Fermi momentum kp for each neutron spin is calculated from 
the density of neutrons in the periodic cube, 

kp = y{3'K^NY^\ (99) 
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FIG. 5: Plots of the three energy ratios defined in Eq. (j95p versus EucUdean projection time t. 
These are labeUed as LO3, ANLO3, NLO3 respectively. 



X. DISCUSSION 

A. Comparisons w^ith other results 

In Fig. [6] we compare the ground state energy ratio Eq/Eq'^'^ for LO2, NLO2, LO3, and 
NLO3 as a function of kp- We note two points here. First the difference between LO3 and 
NLO3 values for Eq/ Eq'^'^ is relatively small over the range of kp plotted. This suggests that 
the convergence of the effective field theory expansion appears reliable, and the difference 
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TABLE IV: Fit results for Sq.lOs /^o'^'' and ^cNLOg/^j 



1 V 




t\j p I iVlc V ) 


TP 1 TTifree 
-C^O.LOs / -f^O 


T? 1 Tplree 
-C/O.NLOs/ -C/Q 




O 


A 




4fiQ('9'l 


n 43fil'9'> 


6 


O 




1 94 


^1Q('4'l 




8 


o 


U 




\} .0'0'-±y-± \ 


U.OO / I "4: 1 


u.o 


o 
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oo 




D ^71 ('S'l 




12 


5 


142 
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FIG. 6: Comparison of the ground state energy ratio Eq/Eq'^'^ for LO2, NLO2, LO3, and NLO3 as 
a function of fci?. 



between LO3 and NLO3 values provides an upper estimate on the size of contributions at 
higher orders. Second the results for NLO2 and NLO3 agree for kp less than 100 MeV. 
This is the region where we expect the perturbative treatment of NLO2 corrections to be 
accurate. The agreement with NLO3 provides some confidence in the effective field theory 
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4l[, GC 2007 



free 



421, and GIFPS 2008 



43]. 



for LO3 and NLO3 versus Fermi momentum kp. For 



APR 1998 



391], CMPR v6 and v8' 2003 



40|], SP 2005 



approach to dilute neutron matter. It is also an explicit test of model independence at fixed 
lattice spacing as suggested in Ref. 

In Fig. [7] we compare ground state e nerg ies for LO3 and NLO3 with other results from 



the literature: FP 1981 



APR 1998 



39], CMPR v6 and v8' |40|, SP 2005 |4l|, GC 2007 



42I ]. and GIFPS 2008 j43|. Compared with other calculations our ground state energies are 



slightly lower for kp near 130 MeV, but overall there is relatively good agreement. 



B. Expansion near the unitarity limit 

The unitarity limit is an idealized limit of attractive two-component fermions where the 
S-wave scattering length is infinite and the range of the interaction is negligible. The S- 
wave scattering length for neutron-neutron scattering is —18.5 fm, while the range of the 
interaction is comparable to the Compton wavelength of the pion, m~^. The unitarity limit 
is approximately realized in neutron matter when the average particle separation is between 
these two length scales. This occurs at a Fermi momentum of about 80 MeV. In the 
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unitarity limit the ground state has no dimensionful parameters other than particle density. 
Therefore the ground state energy in the unitarity limit should obey a simple and universal 
relation Eq = C.Eq'^'^ for some dimensionless constant ^. 

The unitarity limit has been reproduced in trapped cold atom experiments using ^Li and 
^°K. The scattering length is tuned to infinity using a Feshbach resonance and the system 
is sufficiently dilute that the range of the interaction is negligible. Recent experimental 
measurements for ^ give 0.32+J° ^ , 0.51(4) OAGtH M, and 0.39(2) llTj. There nave 



been numerous analytic calculations for ^ varying over the range from 0.2 to 0.6 [48 



49 



51 



52 



53 



54 



55 



56 



57 



50, 



581]. Several numerical calculations both on the lattice and in the 



continuum find results varying from about 0.25 to 0.45 



68 



59 



60 



62 



63 



64 



65 



66 



69 



7C 



67|, 



7l| . The most recent of these numerical calculations agree on a smaller window 



between 0.30 and 0.40. 

For finite S'-wave scattering length the deviation away from unitarity can be parame- 
terized as 

As shorthand notation we define 

6 



fikpao) = ^ 



(101) 
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iterature there is gen eral agreement on the value of ^i, ranging from about 



60 



66 



68 



72 



73l |. In the following analysis we use the values ^ = 0.31 



In the recent 
0.8 to 1.0 

and .^1 = 0.81 calculated in Ref. [6. 

In addition to the corrections at finite scattering length we expect corrections proportional 
to kpTf) due to the S'-wave effective range tq. For neutron- neutron scattering tq is 2.7 fm. 
We also expect higher-order corrections away from the unitarity limit arising from higher 
powers of l/^kpao) and kpr^, as well as other terms associated with the S'-wave shape 
parameter and triplet P-wave scattering volumes. In general we can write 



(102) 



Due to the relatively narrow spread oi kp values considered in our lattice simulations, it is 
difficult to constrain C2 and C3 and higher coefficient powers. However we can constrain the 
parameter ci. 
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FIG. 8: Comparison oi Eq^^i^o.J E^^'^ with various fits involving subsets of the unknown parameters 
ci,C2,C3 as defined in Eq. (jl02p . 



If we set C2 = C3 = and determine Ci from the data point with the smallest value for kp 
we get 

^o,NL03/^o"' ~ fikpao) + 0.14fc^ro. (103) 
If instead we set C3 = and determine ci and C2 simultaneously we find 

^CNLOa/^S"' ~ fikpao) + 0.27A;^ro - Q.Uklm-\ (104) 

As a third alternative if we set C2 = and fit ci and C3 simultaneously, we get 

E^^^i^ojE'r - fikpao) + O.lTfc^ro - 0.264m;l (105) 

The results of these fits are shown in Fig. [HI Our simple analysis suggests a value for ci 
in the range between 0.14 and 0.27. This is consistent with the value 0.15 for the same 
coefficient found in Ref. p]. 

XI. SUMMARY 

We have presented lattice simulations for the ground state energy of dilute neutron matter 
at next-to-leading order in chiral effective field theory. We have solved some problems that 
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arose in recent work using leading-order lattice actions LOi and LO2. LOi involved point- 
like "contact" interactions while LO2 used Gaussian-smeared "contact" interactions. In this 
work we introduced a new action LO3 which equals LO2 in each 5'-wave channel and equals 
LOi in each P-wave channel. The action was constructed using projection operators for 
the spin-singlet/isospin-triplet and spin-triplet/isospin-singlet channels. Using the spherical 
wall method we computed phase shifts and mixing angles for the new lattice action up to 
next-to-leading order and fitted all unknown operator coefficients. 

In the auxiliary-field formalism we used Euclidean-time projection Monte Carlo to com- 
pute the ground state energy oi N = 8, 12, 16 neutrons in a periodic cube, covering a 
density range from 2% to 10% of normal nuclear density. For kp less than 100 MeV we 
found ground state energies at next-to-leading order that agreed with earlier lattice results 
using the action NLO2. For kp greater than 100 MeV we found that the new action leads 
to much smaller corrections at next-to-leading order. The difference between leading-order 
and next-to-leading-order values provides an upper estimate on the size of contributions at 
higher orders. Though we find somewhat lower values for the ground state energy near 
kp = 130 MeV, our results are in general agreement with other calculations reported in the 
literature. 

The ground state energy ratio Eq/Eq^'^ was also analyzed as an expansion about the 
unitarity limit. We considered corrections due to finite scattering length, nonzero effective 
range, and higher-order effects. If we use the parameterization 

Eq/E^^^" ^ fikpao) + cikpro + C2kpmf + c^klm,-^ + • • • , (106) 

we find ci in the range from 0.14 to 0.27. In principle the coefficient ci is a universal 
constant that can be measured in any two-component fermionic system near the unitarity 
limit. Explicit tests of this universality may be a subject for future investigation. With 
regard to further investigations of neutron matter, future work on lattice simulations should 
be considered at next-to-next-to-leading order in chiral effective field theory. Simulations 
should also be done at smaller and larger lattice spacings to check independence on the 
lattice spacing and to probe both higher and lower densities. 
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